home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C / Utilities / Calc 1.24.7 / lib / surd.cal < prev    next >
Encoding:
Text File  |  1992-02-24  |  4.9 KB  |  297 lines  |  [TEXT/R*ch]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Calculate using quadratic surds of the form: a + b * sqrt(D).
  7.  */
  8.  
  9. obj surd {a, b};        /* definition of the surd object */
  10.  
  11. global surd_type;        /* type of surd (value of D) */
  12. global surd__;            /* example surd for testing against */
  13.  
  14. surd_type = -1;            /* default */
  15. obj surd surd__;        /* set object */
  16.  
  17.  
  18. define surd(a,b)
  19. {
  20.     local x;
  21.  
  22.     obj surd x;
  23.     x.a = a;
  24.     x.b = b;
  25.     return x;
  26. }
  27.  
  28.  
  29. define surd_print(a)
  30. {
  31.     print "surd(" : a.a : ", " : a.b : ")" :;
  32. }
  33.  
  34.  
  35. define surd_conj(a)
  36. {
  37.     local    x;
  38.  
  39.     obj surd x;
  40.     x.a = a.a;
  41.     x.b = -a.b;
  42.     return x;
  43. }
  44.  
  45.  
  46. define surd_norm(a)
  47. {
  48.     return a.a^2 + abs(surd_type) * a.b^2;
  49. }
  50.  
  51.  
  52. define surd_value(a, xepsilon)
  53. {
  54.     local    epsilon;
  55.  
  56.     epsilon = xepsilon;
  57.     if (isnull(epsilon))
  58.         epsilon = epsilon();
  59.     return a.a + a.b * sqrt(surd_type, epsilon);
  60. }
  61.  
  62. define surd_add(a, b)
  63. {
  64.     local x;
  65.  
  66.     obj surd x;
  67.     if (!istype(b, x)) {
  68.         x.a = a.a + b;
  69.         x.b = a.b;
  70.         return x;
  71.     }
  72.     if (!istype(a, x)) {
  73.         x.a = a + b.a;
  74.         x.b = b.b;
  75.         return x;
  76.     }
  77.     x.a = a.a + b.a;
  78.     x.b = a.b + b.b;
  79.     if (x.b)
  80.         return x;
  81.     return x.a;
  82. }
  83.  
  84.  
  85. define surd_sub(a, b)
  86. {
  87.     local x;
  88.  
  89.     obj surd x;
  90.     if (!istype(b, x)) {
  91.         x.a = a.a - b;
  92.         x.b = a.b;
  93.         return x;
  94.     }
  95.     if (!istype(a, x)) {
  96.         x.a = a - b.a;
  97.         x.b = -b.b;
  98.         return x;
  99.     }
  100.     x.a = a.a - b.a;
  101.     x.b = a.b - b.b;
  102.     if (x.b)
  103.         return x;
  104.     return x.a;
  105. }
  106.  
  107.  
  108. define surd_inc(a)
  109. {
  110.     local    x;
  111.  
  112.     x = a;
  113.     x.a++;
  114.     return x;
  115. }
  116.  
  117.  
  118. define surd_dec(a)
  119. {
  120.     local    x;
  121.  
  122.     x = a;
  123.     x.a--;
  124.     return x;
  125. }
  126.  
  127.  
  128. define surd_neg(a)
  129. {
  130.     local    x;
  131.  
  132.     obj surd x;
  133.     x.a = -a.a;
  134.     x.b = -a.b;
  135.     return x;
  136. }
  137.  
  138.  
  139. define surd_mul(a, b)
  140. {
  141.     local x;
  142.  
  143.     obj surd x;
  144.     if (!istype(b, x)) {
  145.         x.a = a.a * b;
  146.         x.b = a.b * b;
  147.     } else if (!istype(a, x)) {
  148.         x.a = b.a * a;
  149.         x.b = b.b * a;
  150.     } else {
  151.         x.a = a.a * b.a + surd_type * a.b * b.b;
  152.         x.b = a.a * b.b + a.b * b.a;
  153.     }
  154.     if (x.b)
  155.         return x;
  156.     return x.a;
  157. }
  158.  
  159.  
  160. define surd_square(a)
  161. {
  162.     local x;
  163.  
  164.     obj surd x;
  165.     x.a = a.a^2 + a.b^2 * surd_type;
  166.     x.b = a.a * a.b * 2;
  167.     if (x.b)
  168.         return x;
  169.     return x.a;
  170. }
  171.  
  172.  
  173. define surd_scale(a, b)
  174. {
  175.     local    x;
  176.  
  177.     obj surd x;
  178.     x.a = scale(a.a, b);
  179.     x.b = scale(a.b, b);
  180.     return x;
  181. }
  182.  
  183.  
  184. define surd_shift(a, b)
  185. {
  186.     local    x;
  187.  
  188.     obj surd x;
  189.     x.a = a.a << b;
  190.     x.b = a.b << b;
  191.     if (x.b)
  192.         return x;
  193.     return x.a;
  194. }
  195.  
  196.  
  197. define surd_div(a, b)
  198. {
  199.     local x, y;
  200.  
  201.     if ((a == 0) && b)
  202.         return 0;
  203.     obj surd x;
  204.     if (!istype(b, x)) {
  205.         x.a = a.a / b;
  206.         x.b = a.b / b;
  207.         return x;
  208.     }
  209.     y = b;
  210.     y.b = -b.b;
  211.     return (a * y) / (b.a^2 - surd_type * b.b^2);
  212. }
  213.  
  214.  
  215. define surd_inv(a)
  216. {
  217.     return 1 / a;
  218. }
  219.  
  220.  
  221. define surd_sgn(a)
  222. {
  223.     if (surd_type < 0)
  224.         quit "Taking sign of complex surd";
  225.     if (a.a == 0)
  226.         return sgn(a.b);
  227.     if (a.b == 0)
  228.         return sgn(a.a);
  229.     if ((a.a > 0) && (a.b > 0))
  230.         return 1;
  231.     if ((a.a < 0) && (a.b < 0))
  232.         return -1;
  233.     return sgn(a.a^2 - a.b^2 * surd_type) * sgn(a.a);
  234. }
  235.  
  236.  
  237. define surd_cmp(a, b)
  238. {
  239.     if (!istype(a, surd__))
  240.         return ((b.b != 0) || (a != b.a));
  241.     if (!istype(b, surd__))
  242.         return ((a.b != 0) || (b != a.a));
  243.     return ((a.a != b.a) || (a.b != b.b));
  244. }
  245.  
  246.  
  247. define surd_rel(a, b)
  248. {
  249.     local x, y;
  250.  
  251.     if (surd_type < 0)
  252.         quit "Relative comparison of complex surds";
  253.     if (!istype(a, surd__)) {
  254.         x = a - b.a;
  255.         y = -b.b;
  256.     } else if (!istype(b, surd__)) {
  257.         x = a.a - b;
  258.         y = a.b;
  259.     } else {
  260.         x = a.a - b.a;
  261.         y = a.b - b.b;
  262.     }
  263.     if (y == 0)
  264.         return sgn(x);
  265.     if (x == 0)
  266.         return sgn(y);
  267.     if ((x < 0) && (y < 0))
  268.         return -1;
  269.     if ((x > 0) && (y > 0))
  270.         return 1;
  271.     return sgn(x^2 - y^2 * surd_type) * sgn(x);
  272. }
  273.  
  274. global lib_debug;
  275. if (!isnum(lib_debug) || lib_debug>0) print "obj surd {a, b} defined"
  276. if (!isnum(lib_debug) || lib_debug>0) print "surd(a, b) defined"
  277. if (!isnum(lib_debug) || lib_debug>0) print "surd_print(a) defined"
  278. if (!isnum(lib_debug) || lib_debug>0) print "surd_conj(a) defined"
  279. if (!isnum(lib_debug) || lib_debug>0) print "surd_norm(a) defined"
  280. if (!isnum(lib_debug) || lib_debug>0) print "surd_value(a, xepsilon) defined"
  281. if (!isnum(lib_debug) || lib_debug>0) print "surd_add(a, b) defined"
  282. if (!isnum(lib_debug) || lib_debug>0) print "surd_sub(a, b) defined"
  283. if (!isnum(lib_debug) || lib_debug>0) print "surd_inc(a) defined"
  284. if (!isnum(lib_debug) || lib_debug>0) print "surd_dec(a) defined"
  285. if (!isnum(lib_debug) || lib_debug>0) print "surd_neg(a) defined"
  286. if (!isnum(lib_debug) || lib_debug>0) print "surd_mul(a, b) defined"
  287. if (!isnum(lib_debug) || lib_debug>0) print "surd_square(a) defined"
  288. if (!isnum(lib_debug) || lib_debug>0) print "surd_scale(a, b) defined"
  289. if (!isnum(lib_debug) || lib_debug>0) print "surd_shift(a, b) defined"
  290. if (!isnum(lib_debug) || lib_debug>0) print "surd_div(a, b) defined"
  291. if (!isnum(lib_debug) || lib_debug>0) print "surd_inv(a) defined"
  292. if (!isnum(lib_debug) || lib_debug>0) print "surd_sgn(a) defined"
  293. if (!isnum(lib_debug) || lib_debug>0) print "surd_cmp(a, b) defined"
  294. if (!isnum(lib_debug) || lib_debug>0) print "surd_rel(a, b) defined"
  295. if (!isnum(lib_debug) || lib_debug>0) print "surd_type defined"
  296. if (!isnum(lib_debug) || lib_debug>0) print "set surd_type as needed"
  297.